Density-functional theory calculations for 
poly-atomic systems: Electronic structure, 
static and elastic properties and ab initio 
molecular dynamics 



Michel Bockstedte, ^Alexander Kley, Jorg Neugebauer and Matthias Schefher 



0^ 
0^ . 

Fritz- Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 
D- 14195 Berlin- Dahlem, Germany 

^ Abstract 
<^ 

I \ The package fhi96md is an efficient code to perform density- functional theory total- 
^ ■ energy calculations for materials ranging from insulators to transition metals. The 
package employs first-principles pseudopotentials, and a plane-wave basis-set. For 
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fSJ I exchange and correlation both the local density and generalized gradient approxima- 
in tions are implemented. The code has a low storage demand and performs efficiently 
on low budget personal computers as well as high performance computers. 
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PROGRAM SUMMARY 

Title of the Program: fhi96md 
catalogue number: 

Program obtainable from: CPC Program Li- 
brary, Queen's University of Belfast, N. Ireland 

Licensing provisions: Persons requesting the 
program must sign the standard CPC non- 
profit use hcense (see license agreement printed 
in every issue). 

Computer for which the program is designed 
and others on which it has been tested: IBM 
RS/6000, Pentium-PC 

Installation: Fritz-Haber-Institut der Max- 
Planck-Gesellschaft, Berlin-Dahlem 

Operating system: UNIX 

Programming language: FORTRAN 77 

Memory required to execute with typical data: 
1-32 Mwords 

No. of bits in a word: 32 

Memory required for test run: 0.5 MB 

Time for test run: 6 min 

No. of lines in distributed program, including 
test data, etc.: 14000 

Keywords: density-functional, local-density, 
generalized gradient, pseudopotentials, plane- 
wave basis, super cell, molecular dynamics, struc- 
ture optimization, total-energy, potential-energy 
surfaces, chemical binding, diffusion, surface 
reactions, crystals, defects in crystals, surfaces, 
molecules 



Nature of the physical problem 

In poly-atomic systems as for example molecu- 
les [1,2], crystals [3,4], defects in crystals [5,6], 
surfaces [7-9], it is highly desirable to perform 
accurate electronic structure calculations, with- 
out introducing uncontrollable approximations. 
This enables the predictive description of equi- 
librium properties as well as of non-equilibrium 



phenomena for a wide class of materials. Ex- 
amples studied with the present code or its 
predecessor include mcta-stabilities of defects 
[5,10], surface reconstructions [11,12], diffusion 
[9,13], surface reactions [14-16], and phase tran- 
sitions [17]. Molecular dynamics simulations 
combined with first- principles forces are a pow- 
erful tool to analyze the motion of the nu- 
clei [13,18] and to accurately calculate thermo- 
dynamic properties such as diffusion constants 
and free energies [13]. 

The computer code described below enables 

this variety of investigations. It employs density- 
functional theory [20] together with the local- 
density approximation [21,22] or generalized 
gradient approximations [23-25] for the exchange- 
correlation functional. 

Method of solution 

Ab initio molecular dynamics on the Born-Oppen- 
heimcr surface is implemented by a two step 
method. In the first step the Kohn-Sham equa- 
tion [26] is solved self-consistently to obtain 
the electron ground state and the forces on the 
nuclei. In a second step these forces are used to 
integrate the equations of motion for the next 
time step. The calculation of the total-energy 
and the Kohn-Sham operator in a plane-wave 
basis-set is done by the momentum-space method [27] . 
To solve the Kohn-Sham equation the pack- 
age fhi96md employs the iterative schemes of 
Williams and Soler [28] and Payne et al. [29]. 
We use the frozen-core approximation and re- 
place the effect of the core electrons by norm- 
conserving pseudopotentials [30-33] in the fully 
separable form [34]. The equations of motion of 
the nuclei are integrated using standard schemes 
in molecular dynamics such as the Verlet algo- 
rithm. Optionally an efficient structure opti- 
mization can be performed by a second order 
algorithm with a damping term. 

Restrictions on the complexity of the problem 

Only pseudopotentials with s-, p-, and (i-compo- 
nents are implemented. For highly correlated 
systems (e.g. /-electrons) the treatment of the 
exchange-correlation interaction is not appro- 
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priate. Relativistic effects are included via scalar 
relativistic pseudopotentials. The system is as- 
sumed to be non- magnetic, but a generaliza- 
tion of the program to magnetic states is straight 
forward. 
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LONG WRITE-UP 



1 Introduction 

Total-energy calculations and molecular dynamics simulations employing density-functional the- 
ory [20] represent a reliable tool in condensed matter physics, material science, chemical physics 
and physical chemistry. A large variety of applications in systems as different as molecules [1,2], 
bulk materials [3-6] and surfaces [7-9] have proven the power of these methods in analyzing as 
well as in predicting equilibrium and non-equilibrium properties. Ab initio molecular dynamics 
simulations enable the analysis of the atomic motion and allow the accurate calculation of ther- 
modynamic properties such as the free energy, diffusion constants and melting temperatures of 
materials. 

The package fhi96md described in this paper is especially designed to investigate the material 
properties of large systems. It is based on an iterative approach to obtain the electron ground 
state. Norm-conserving pseudopotentials [30-33] in the fully separable form of Kleinman and 
Bylander [34] are used to describe the potential of the nuclei and core electrons acting on the 
valence electrons. Exchange and correlation are described by either the local-density approxima- 
tion [21,22] or generalized gradient approximation [23-25]. The equations of motion of the nuclei 
are integrated using standard schemes in molecular dynamics. Optionally an efficient structure 
optimization can be performed by a damped dynamics scheme. 

The package fhi96md is based on a previous version fhi93cp [35]. Advances of the new version 
are an improved iteration scheme to solve the Kohn-Sham equations, the generalized gradient ap- 
proximation for exchange and correlation, a mixed basis-set initialization and molecular dynamics. 
The package consists of the program fhi96md and a start utility start. The program fhi96md can 
be used to perform static total energy calculation or ab initio molecular dynamics simulation. The 
utility start assists in generating the parameter file and the input file required to compile and run 
fhi96md, thereby utilizing a low memory demand for each individual run. 

The paper is organized as follows. In Section 2 we bricfiy outline the method as implemented 
in the package. Section 3 describes the program structure and the input and output files of the 
program. The last two Sections concern the installation of the package and the test run. Input 
and output of the test run are found at the end of the paper. Appendix A hsts the formulas and 
expressions implemented in the package and Appendix B describes the parameter and input files 
as generated by the start utility. 



2 Ab initio molecular dynamics 

The key approach to molecular dynamics (MD) is (i) the separation of the slow motion of the 
nuclei from the fast motion of the electrons within the Born-Oppenheimer approximation and (ii) 
that the motion of the nuclei is governed by Newton's equation of motion: 




dRj 



d 



Eo ({R.}) 



(1) 
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where Mj and Rj are the masses and coordinates of the N^^i atoms rsp. and Eq ({Rj}) is the 
many-electron ground state energy. We employ density-functional theory (DFT) together with 
common approximations to the exchange-correlation energy functional to obtain the ground state 
electron density and the forces acting on the nuclei. Computationally the ground state electron 
density is obtained by self-consistently solving the Kohn-Sham equations in a pseudopotential 
plane-wave approach. Once the forces acting on the nuclei are calculated the equations of motion 
are integrated numerically. To perform simulations in the canonical ensemble a Nose-Hoover 
thermostat [36,37] is employed to control the temperature during the simulation. The methods 
involved in finding the electronic ground state and integrating the equation of motion of the nuclei 
are outlined in the first two Subsections. 

Frequently one is only interested in the geometry of stable and metastable structures, e.g. of 
surfaces, defects, or complicated crystals. The equilibrium geometry and the corresponding total- 
energy can be obtained by relaxing the coordinates of the nuclei starting from a guessed geometry. 
An example for this class of application is the calculation of adiabatic potential energy surfaces, 
which are widely used to study surface reactions or defect migration on the microscopic level. An 
outline of the techniques for structure optimization is given at the end of this Section. 

2.1 Solving the Kohn-Sham equations 

2.1.0.1 The energy functional The key variable in DFT is the electron density n(r). As 
stated by the fundamental theorem of Hohenberg and Kohn [20] the ground state energy Eo{{Iij}) 
of the system for given positions of the nuclei {Rj} is the minimum of the Kohn-Sham total energy 
functional [26] with respect to the electron density n. The total energy functional E [n] is: 

E [n] = [n] + E^ [n] + E""" [n] + E^^ [n] + , (2) 

where is the kinetic energy of non-interacting electrons, E^ is the Hartree energy, and E'^'^ 
is the exchange-correlation energy. The energy expressions are briefly discussed in the following. 
The explicit expressions for each of the contributions to the total-energy, potentials and forces as 
implemented in the program can be found in Appendix A. 

The energy of the electron-nuclei and nuclei-nuclei interaction and £;°^c-nuc 

£;--"^[n] = /(iVV^"-""'=(r)n(r) and £;nuc-nuc ^ 1 ^ , (3) 

where Zj and Zj are the charges of the corresponding nuclei. Throughout the paper we employ 
atomic units (energy in hartree) unless otherwise noted. As approximations to the exchange- 
correlation energy functional E^*^ [n] we employ the local-density approximation (LDA) - as 
obtained from the homogeneous electron gas by Ccplerley and Alder [21] in the parameterization of 
Perdew and Zunger [22] - and the generalized gradient approximations of Becke and Perdew [23,24] 
(BP), and of Perdew et al. [25] (PW91). 

The system is represented by the super cell approach. The geometry of the nuclei is contained in 
a super cell, which is periodically repeated on a lattice. The coordinates Rj of a nucleus or its 
periodic image hence are 

Rj^T/sA + R- with J = {4,4,R} , 
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where the index 4 for the species of a nucleus, the index /a for the nucleus itself and the lattice 
vector R pointing to the origin of the cell or of its image. 

The effect of the core electrons and the coulomb potentials of the nuclei is replaced by soft 
pseudopotentials which enables the efficient use of a plane-wave basis: 

ye-nuc ^ ^ ^ ^^e-nuc^^ _ ^^^^^^ _ R) (4) 

Wc employ norm conserving pseudopotentials constructed e.g. following the schemes of Hamann [31] 
or Troulher and Martins [32] . These have proven to yield transferable potentials for a broad class of 
nuclei ranging from first row elements to transition metals. The pseudopotentials are represented 
in the fully separable form as proposed by Kleinman and Bylander [34] 

imax I |AT/nl|/PS \/./,PS |A\/"M 
^h,hoc^ 2^ /^/,ps |AT/nlU/,ps \ ' 



1=0, rn=-l 



{'^^s,l,m\^^ll]l\'^h,l,m) 



where AVj^J;(r) = V/^,i(r) — Vj^^i^^^{r), and Vj^^i{r) are the radial components of the semi- local 
pseudopotential and ^^^^^^^(r) = R\''{r)Yi^{9, 0) are the node-free atomic pseudo wave functions. 
In this form the pseudopotential is sphtted into a local part V}P^'^°^^' and a non-local but separable 
part V}P^'°'. Correspondingly the potential V^~™^^ and the energy E^"™^^ are expressed as 

ye-nuc _ yps,Iocal _|_ ^ps,nl ^e-nuc _ ^ps,local _|_ ^ps,nl ^g-j 



For some atomic species such as sodium, zinc and copper, it is essential to include the non-linear 
core valence exchange- correlation description (NLCV-XC) [38,39]. In this case the effect of the 
core electron density on the exchange-correlation energy is described by a pseudo core density 
ff°^'^{Y), which is the superposition of the smoothed core densities hfJ^{\Y — Ti^j^ — 'R\) constructed 
together with the pseudopotentials [40]. The exchange- correlation functionals E^'^ [n] and [n\ 
are then replaced by E^^ [n + n^'^^^] and [n + ff^^^] . 



2.1.0.2 Kohn-Sham equations In the Kohn-Sham scheme [26] the electron density is ex- 
pressed by a set of orthogonal, normalized Kohn-Sham orbitals 0Q;(r) 

n{v)^Y.fMv)\\ (7) 

a 

The occupation numbers vary between and 2 as the electron spin is not included explicitly 
and the sum over all occupation numbers is the total number of electrons N^i per super cell. The 
ground state electron density is calculated by solving the Kohn-Sham equations self-consistently 
for these orbitals 

(-^V^ + T/^-- + + V"^' [n]) 0,(r) = 6« 0,(r) . (8) 

HKS 

This is equivalent to a constrained minimization of the total energy functional E[n\ with respect 
to the Kohn-Sham orbitals 

£"0 = min E[n] with {4'a\4'i3) = / '^(r) d'^r = N^i- (9) 

{\<t>a}} J 
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2.1.0.3 Occupation numbers The occupation numbers are determined by the Fermi 
distribution 

fa = 2^^ , (10) 

e '=s^ci + 1 

where the Fermi energy /i/ is defined by the number of electrons N^i = J2afa and Td is the 
temperature. It is convenient to use an artificially high temperature (e.g. k-sTgi = 0.1 eV), since 
a broadening of the occupation improves the stability and the speed of convergence of the cal- 
culation. However, as pointed out in Ref. [41] this corresponds to minimizing the free energy 
functional F[n] = E[n] — Td with the entropy 5*61 — J2aifa^^fa + (1 ~ fa)^^{i — fa)) rather 
than the total-energy E[n]. In order to obtain the total-energy at Td = OK, we employ that 
£"0 = E[n] — |Td S'd yields the total-energy at Td = OK up to 0(T^i) [8,42]. However, for typical 
values of ksT^i the deviations are small and the forces on the nuclei remain accurate [43]. 



2.1.0.4 Plane-wave basis-set The Kohn-Sham orbitals are represented by a plane-wave 
basis-set 

0.k(r)= E Q,G+ke^(«+'^)-- (11) 

G,i|G+k|2<Ecut 

truncated at a energy cut-off £^cut- The Brillouin zone integral over the k-points is replaced by 
a sum over a set of special k-points with the corresponding weights [44]. In the start utility 
the Monkhorst-Pack scheme [44] has been implemented to generate special k-points and for each 
k-point set the quality check proposed by Chadi and Cohen [45] is automatically performed. 



2.1.0.5 Solving the Kohn-Sham equation In the past few years iterative techniques have 
become the method of choice to solve the Kohn-Sham equations and enabled first-principles 
studies even for large systems. The key idea is to minimize the energy functional with respect 
to the wave function |0i,k) starting with a trial wave function The energy minimization 

scheme is formulated in terms of an equation of motion (EOM) for the wave function The 
simplest scheme is the steepest descent approach 

^i4j>=(ev-^Ks)i0a> ' (12) 

under the ortho-normality constraint (0i^|0j^k) = where the e^^k are Lagrange parameters 
due to the ortho-normality constraint. In the program a more refined and efficient scheme based 
on a second order EOM has been implemented 

^I0S> + 27 = (e.,k - ^Ks) , (13) 

where 7 is a damping parameter. The EOM is integrated for a step length St by the Joannopoulos 
approach [29], which iteratively improves the initial wave functions. Tassone et al. [47] noted that 
such a scheme is as efficient as the conjugate gradient techniques [48]. Though conjugate gradient 
techniques need a smaller number of iterations to minimize the energy functional, they need 
to perform an accurate line minimization and have to calculate the conjugate search directions 
making this algorithm rather costly compared to the calculation of the steepest descent direction. 
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Hence, the effort to perform an iteration of a conjugate gradient algorithm is substantially larger 
than to perform an iterative step using the damped Joannopoulos algorithm. 

In the damped Joannopoulos algorithm the new wave function is constructed from the 

wave functions of the last two iteration steps r and r — 1 



(G + k|0l;-+^)) = (G + k|0g) +/?G (G + kl^S) (14) 

+ 7G (G + kl^;;^"'^) - (G + k|^Ks|0S5) 



where the coefficients are 



G 

7G = e 
Vg 



ei,y,{h{St) - 1) ^ (G 4 


-k|iJKs 


G + k) e-T"^* 


ei,k- (G + k|^KS 


G + k) 



(15) 
(16) 



h{St) - e-^^^ - 1 



Q,k-(G + k|i/Ks|G + k) 

with the damping parameter 7, the step length parameter 5t, and with Cj ^ 
function h{5t) is defined by 



;^I^Ks|0S2).The 



h{dt) 



e 2 cos (u 5t) 



if a;2 > 



e-i-^* cosh {^/\^\St) ifu^ <0 



with = ei,k - (G + k|^Ks|0tk) - T- 

After each iteration step the wave functions {|0i,k)} have to be ortho-normalized, which is done 
by the Grahm-Schmidt ortho-normalization scheme. Otherwise all states would converge to the 
lowest lying state. 

The steepest descent direction which points to the total-energy minimum is given by 

(G + k|^^S|^,,k) = (G + k| - ^V>,,k) (17) 

+ (G + k| t>H + yP^,l°cal ^ yXC 1^.^^^ 

y'local 

+(G + k|yP^'-Vi,k) . 

It is calculated very efficiently avoiding expensive vector matrix products by evaluating 
(G + k| — |V^|0j,k) in the momentum representation and (G + k|y^°^^'|0j k) in configuration 
space representation where and are diagonal. Transformations between momentum 

representation and configuration space representation are performed by fast Fourier transforma- 
tions and the cost for the calculation of the latter two expressions is 0{N In N) operations per 
state and k-point [46]. 
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Though the damped Joannopoulos algorithm is at least twice as efficient as the first order scheme, 
given by Eq. (12), additional storage for the wave function 10-^"^^) is needed. Hence we use the 
William-Soler algorithm [28] whenever storage requirements do not permit to employ the damped 
Joannopoulos algorithm. The coefficients of this scheme in Eq. (14) are 

~ 6i,k-(G+k|HKs|G+k> 

e,,k-(G + k|^Ks|G + k) 

(■>r,,k-(G+k'i}KH'G+k\ 

''''~ei,k-(G + k|^Ks|G + k) 

with 7g = 0. The wave functions are improved by successive iterations by either of the schemes 
outlined above until a convergence criterion concerning the accuracy of the energy or the forces 
as described in Section 3.1 is fulfilled. On convergence the package proceeds with a MD step, a 
structure optimization step or just terminates. 



2.1.0.6 Mixed-basis-set initialization As pointed out above iterative techniques require 
an initial guess for the wave function 10^°^). Regardless of the method used, a good initial guess 

for the initial wave function is essential and can significantly improve the convergence of the 

method. The simplest choice is to generate 10^'^) from random numbers. However, the trial wave 
function obtained by this procedure will be far away from the final solution. Consequently, less 
iterations are necessary when the initial wave function is obtained in a more physical way, i.e. by 
a direct diagonalization of the Kohn-Sham Hamiltonian. The effort to diagonalize the Kohn-Sham 
Hamiltonian on the full basis, however, outweighs the effort saved in the iterative diagonalization. 
Instead, a very efficient procedure is to represent the Kohn-Sham orbitals {|x!^,k)} in a LCAO or 
mixed-basis-set which dramatically reduces the number of basis functions and at the same time 
gives a good approximation to the solution within the full basis [49]. In this basis the Kohn-Sham 
Hamiltonian is easily diagonalized for a few number of self-consistency cycles and the overall effort 
is considerably reduced. 

The mixed basis-set is a subset of the full plane-wave basis. It includes Bloch-states derived 
from atomic orbitals and plane waves up to an energy cut-off E"^^ much lower than the energy 
cut-off of the full basis. The locahzed basis functions are constructed by Grahm- Schmidt ortho- 
normalization from the states lili^j^^i^m) , which are projections of Bloch-sums onto the plane- wave 
basis: 

\fiiMm)= E (G + k|/.,,,,,,,,J|G + k) . (18) 

G 

£i^\f<^lG+k|2<iS<;ut 

The functions iJ,i^j^,i,m{r) are 

l^iMmir) = Ee''-''C^,m(r - R - T,,,/J , (19) 

R 

where ip'i^imi^) atomic pseudo wave functions. Using this basis an adequate description 

of the localized states as well as of the interstitial region is possible and the adequate description 
of strongly localized states reduces the number of necessary iterations of the iterative energy 
minimization scheme [49]. 
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In Eq. (18) the localized basis functions are defined in terms of a plane- wave representation. 
Therefore the mixed basis matrix representation {Xti,k\HKs\Xu,k) of the Kohn-Sham operator can 
be easily calculated using the same routines as employed to compute the steepest descent vector 
in Eq. (17) just at the additional cost of a scalar product: 

{Xf.,k\HKs\x.M) = E(XM,k|G + k)(G + k\HKs\Xu,k) ■ (20) 

G 

This fact combined with the small size of the basis reduces the computational effort of the matrix 
diagonalization approximately to the cost of an energy minimization step within the full basis. 

This procedure is repeated until electronic self-consistency is reached. Old and new electron 
densities are mixed using the Broyden-scheme [50]. The resulting Kohn-Sham orbitals are used 
as initial wave functions for the iterative energy minimization scheme described above 



2.1.0.7 Implementation The program fhi96md efficiently performs on low budget personal 
computers even for large systems. Obviously, it works even better on high performance computers. 
The high efficiency of the code is achieved by a low storage demand and an economic cache 
utilization. For example single precision arrays are used to store the wave function coefficients 
Ci,G+k and the reciprocal lattice vectors G without effecting the accuracy of the results. The 
portability of the code to a variety of platforms without loosing efficiency is achieved by using 
standard library subroutines like BLAS-routines. In this way on each platform full advantage can 
be taken of optimized, platform specific implementations of the libraries. However, when such 
libraries are not available (e.g. on a PC) we provide the necessary routines in a separate library 
which can be linked to the program. 

In typical applications the calculation of the electron density and the steepest descent direction 
dominate the computational effort (mainly FFT). These contributions scale like 0{N'^liiN) with 
the system size. In very large systems, e.g. a Silicon super cell with 128 atoms and more, the 

calculation of the (non-local contributions of the) atomic forces (which scale like 0{N'^) becomes 
the dominant contribution. The computational effort for this contribution is dramatically reduced 
by calculating the forces only for the last few electronic steps, when the Born-Oppenheimer surface 
is almost reached. 

2.2 Integrating the equations of motion of the nuclei 

Once the ground state of the electrons is calculated as described in Sec. 2.1 the atomic EOMs 

are integrated with standard MD techniques. We have implemented two schemes, the Verlet 
algorithm and a predictor-corrector algorithm. The choice of the scheme depends on the actual 
type of application. For example, the Verlet algorithm is more stable with respect to an energy 
drift than the predictor-corrector algorithm. However, when it is important to obtain very accurate 
velocities, the predictor-corrector algorithm should be employed. 

The Verlet algorithm [51] is 

Ti^,iSt + 5W) = 2Tj^,St) - rj^,iSt - 5Uc) + 5Cc ^^-^-(\7-^-(^^>^ (21) 
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where Mj^ is the mass of the nuclei and Stnuc is the time step. Prom a numerical point of view 
it is desirable to chose a time step Stnuc as large as possible. A good choice for Stnuc is of the 
shortest phonon period in the system, which is approximately one order of magnitude larger than 
the time step used in the Car-Parrinello ab initio MD [46]. The time step there is limited by the 
fast oscillatory motion of the electrons in the fictitious electron dynamics [52], which originates 
from the need to keep the electrons close to the Born-Oppenheimer surface. 

The velocity of a nucleus is calculated by 

= (3T/3,7a(^) - 4T7,,/,(i - Stnuc) + ri^,jSt - 2 Stuuc)) , (22) 

which is correct to second order in St^uc [53]. Note that the velocities are not directly integrated 
by the algorithm, but obtained from the trajectory itself. 

In a predictor- corrector scheme the discretization error is of higher order in St^uc than that of 
the Verlet algorithm. The usual predictor-corrector (PC) schemes employ an extrapolation step 
to predict positions and velocities, which are then corrected in the corrector step. Por small 
time steps the accuracy of the PC scheme is much higher than that of the Verlet algorithm. 
However, for a large time step it looses this advantage [54], since the error in the extrapolation 
step increases strongly with an increasing time step. Therefore we have implemented a predictor- 
corrector scheme which avoids the extrapolation. In this scheme an Adams-Bashforth predictor 
step 



^h,lS^n) = T /„i^(tn-l) + Stuuc ^ Ciq^k'^ i^j^{tn-k) 

k=l 

and (23) 

,.0 \_^, (J. ^ I ^, ^Js,Ja({^Js,Ja(^n-fc)}) 

fe=l "^»s 



is followed by an Adams-Moulton corrector step 



T hjS^n) = "T h,ia.{'tn-l) + St^uc Pq,0 ^%,iS'^n) + St^uc X] Pq,k^ h,iS^n-k) 

k=l 

and (24) 

^h,lA'^n) — V7„i3,(tn-lj + Otnuc Pqfi 

a Fj^.iad'T Js,/a(^n-fc)}) 

with the coefficients aq^k and f3q^k as tabulated in Ref. [55, pp 126] and with t„ = t+n 5tnuc- Though 
in this scheme the forces have to be re-calculated after the corrector step, this is inexpensively 
done as the corrector step corrects the positions of the nuclei only slightly and convergent forces 
are typically obtained within a few iterative steps. Thus the accuracy of the predictor-corrector 
scheme does not suffer from a large time step. Moreover, the velocities are smooth functions of time 
and the total-energy of the coupled electron nuclei system is practically free from fluctuations. 
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Energy losses encountered in long simulations within the micro-canonical ensemble can be com- 
pensated by a periodic rescaling of the atomic velocities according to a preset average kinetic 
energy per particle. Temperature is optionally enabled by a Nose-Hoover thermostat [36,37]. One 
additional degree of freedom is employed to simulate the heat bath and an extra force is accel- 
erating or de-accelerating the nuclei to maintain the temperature of the system. The modified 
EOMs of the extended system are 



and (25) 
Q = E Miyj^j^ - 9 kuT 

where g is the number of independent degrees of freedom in the system. 

The initial coordinates and velocities of all independent degrees of freedom completely determine 
the trajectory. The initial velocities are generated from Gaussian distributed random numbers 
such that the center of mass is at rest and that the kinetic energy of the system corresponds to 
an initial temperature. In a micro-canonical ensemble the initial temperature together with the 
potential energy of the initial configuration of the nuclei determine the resulting average kinetic 
energy per particle. 

Both schemes, the Verlet algorithm and the predictor-corrector scheme, require knowledge of 
the trajectory at previous time steps. Therefore, to start these algorithms, the Runge-Kutta 
scheme [55] is used to integrate the equation of motion for the first few time steps. 

2.3 Structure optimization 

An efficient and numerically stable method to find the equilibrium geometry is damped Newton 
dynamics. Starting from an initial configuration the nuclei are moved according to the iterative 
scheme 

tSD:^^ = (1 + A - A,, rix'' + l^i. F^,^ ({rK}) , (26) 

where A/3 and /ij^ are the damping and reciprocal mass parameters rsp. The parameters Aj^ and 
/x/g determine whether the nuclei loose their initial potential energy slowly in an oscillatory-like 
motion or whether they move straight into the closest local minimum. The program allows also 
to confine the configuration space open to the search. An example is the calculation of adiabatic 
potential energy surfaces, where the ad-atom or a defect is held fixed and all other atoms are 
allowed to relax. 



3 The package fhi96md 

The package contains the program fhi96md and a start utility start. The program fhi96md per- 
forms the MD simulations and the total-energy calculations. The start utihty generates the file 
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para meter, h containing the parameters and the file input.ini containing the input data, which are 
necessary to compile and to run the program fhi96md. Features of the start utility include an 

automatic search for point group symmetries and the symmetry center in the system. Further, it 
automatically optimizes the parameter file to minimize the memory demand for each individual 
run. 

The following two Subsections describe the structure of the program fhi96md and discuss the 
input files processed by the start utility and the output files generated by fhi96md. 



3. 1 The program structure 

The structure of the program fhi96md is sketched in Fig. 3.1. The self-consistent calculation of the 
electron ground state forms the main body of the program, which is displayed on the left-hand 
side of Fig. 3.1. The movement of the atoms is accomplished in the block "move atoms", which is 
sketched on the right hand side of Fig. 3.1. Note, that the generation of output is not exphcitly 
accounted for in the fiowchart and we refer to it at the end of this Subsection. 

The first block in the flowchart is the initialization block, where the program reads the input files 
inp.mod, inp.ini and the pseudopotential data. Then the routines calculating form factors, structure 
factors and phase factors (c.f. Appendix A. 3) are called and the initial wave functions are set up 
either from a restart file or by a few self-consistency cycles using the mixed basis-set initialization 
(c.f. Section 2.1). Having obtained the initial wave function {"^fl) the program enters the self- 
consistency loop. First, the electron density and the contributions to energy, potential and forces 
are calculated. Note, that the forces are only calculated during MD simulations and structure 
optimization when the electrons are sufficiently close to the Born-Oppcnheimer surface. 

Within the block "move atoms" the atomic EOMs are integrated for one time step in a MD 
simulation or a structure optimization is performed, provided the electrons are sufficiently close 
to the Born-Oppenheimer surface, i.e. the forces are converged. The control over the calculation 
of the forces is handled by this block. If the nuclei have been moved, i.e. either the atomic EOMs 
have been integrated for another time step or one structure optimization step has been executed, 
it also recalculates the structure and phase factors and other quantities that explicitly depend on 
the positions of the nuclei. Upon the first call to the routine f iondyn in this block the restart file 
fort. 20 is read, if provided, and all necessary steps are taken to restart or initialize the dynamics. 

The following two blocks update the wave functions using the damped Joannopoulos or the 
William-Soler algorithm and ortho-normalize the wave function by a Grahm-Schmidt ortho- 
normalization. In the last block within the self-consistency loop the occupation numbers arc up- 
dated e.g. for a metalUc system according to a Fermi distribution by a damped pseudo-eigenvalue 
scheme [35,56]. This block enables also an interactive control over the step length St and damping 
parameter 7 of the energy-minimization scheme, the mass parameter /x/^ of the structure opti- 
mization scheme and the remaining numbers of iterations while the program is running. These 
parameters are updated from the files delt (5r), ionJac (yU/J, stopfile (remaining number of elec- 
tronic iterations) and stopprogram (remaining number of structure optimization steps). If these 
files are empty the parameters are not changed. Finally, the convergence criteria are checked. The 
program terminates when convergence is achieved or when the preset number of iterations or the 
allowed CPU-time is exceeded. 
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Output is generated at the last block of the self-consistency loop and by the routines f iondyn, 
fionsc, fermi and vofrho. The routines f iondyn and o_wave generate restart files for MD 
simulations and total-energy calculations. 

In the mixed-basis-set initialization, the self-consistency loop closely follows the organization of 
that discussed above. First of all the initial electron density is obtained either from a superpo- 
sition of contracted atomic pseudo densities or from an electron density of a previous calcula- 
tion (file fort. 72). The local contributions to the potential and the energy are calculated by the 
routine vofrho. The localized orbitals to construct the mixed-basis-set are set up by routine 
project_init. The non-local contributions to the potential and the energy in the localized basis 
set are calculated by the routine nlrhkb_bO. In the second step the Hamiltonian is constructed 
with the help of routine df orce_bO. The Hamiltonian is diagonalized by standard diagonaliza- 
tion routines. The new electron density is calculated (routine rho_psi). Finally the new electron 
density is mixed with the old density by a Broyden mixing (routine broyden). 

3.2 Input and output files 

Two input files are required as input for the start utihty. The file in p. mod contains the control 
parameters for the run. The file start. inp describes the geometry of the super cell, the configuration 
of the nuclei and parameters relevant for the MD simulation or the structure optimization, and the 
calculation of the electron ground state. The parameters and data contained in the files inp. mod 
and start. inp are described in Tables 3.1 and 3.2. Related entries are grouped in one line as shown 
in the hstings of the input files of the test run in Section TEST RUN. No specific format is 
requested other than that implied by the type of the entries as displayed in the second column 
of Tables 3.1 and 3.2. The start utility processes the files inp. mod and start. inp and generates the 
parameter file parameter. h and the input file inp.ini. These files are described in Appendix B. 

In the following we describe the strategy to set up the files inp. mod and start. inp. First of all 

there are the logical parameters tdyn and tfor- we refer to parameters in the file inp. mod, if not 
noted otherwise. With these two parameters one instructs fhi96md to perform a MD simulation, 
a structure optimization, or if both are set to .false, just to calculate the electronic structure 
for the given configuration of the nuclei. 

For a MD simulation the parameters Ldyn, norder, and delt_ion specify the scheme for integrating 
the equation of motion of the nuclei and the time step (c.f. Section 2.2). Additional parameters in 
the file start. inp determine the simulation ensemble {nthm, Q, T_ion and nfi^rescale ), the set up 
of initial positions and velocities {npos and coordwave - restart options included) and the masses 
of the nuclei M/^. The parameters needed for the structure optimization ///^ and A/^ {ionjac and 
ion-damp) are specified in start. inp as well (c.f. Section 2.3). 

The geometry of the super cell, the positions of the nuclei and optionally the velocities (c.f. 
npos) are specified in the file start. inp. The parameter ihrav and pgind determine the lattice type 
and the point group symmetry of the configuration of the nuclei. With pgind = an automatic 
search for the point group symmetries and the symmetry center is performed. For pgind > 1, 
the symmetry center is the origin (0, 0, 0) of the super cell. In MD simulations and structure 
optimization point group symmetries are usually not applicable. For each of the nsp atomic 
species one declares the properties of the pseudopotential (zv, Lmax and Lloc) and the radius of the 
screening charge n^^"^(r) {rgauss, c.f. Appendix A). The positions of the nuclei tauO and optionally 
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also the velocities vauO follow this declaration (c.f. npos). Note, that lines containing the velocity 
of a nuclei immediately follow the line with the corresponding coordinates. Pseudopotential data 
needs to be provided in the files fort. 11, fort. 12, ... for the nsp atomic species. The data in each 
file is expected either in a parameterized form or tabulated on a logarithmic mesh (parameter 
tpsmesh) . 

Having set up the basic configuration data, we now turn to the data needed for the calculation of 
the electron ground state. The parameter ecut specifies the energy cut-off of the plane-wave basis- 
set. The data xk and wkpt declare the nkpt k-points (LkpoinLrel specifies the frame of reference). 
The start utility reduces the k-point set according to the point group symmetries as specified 
by the parameter pgind. A special k-point set according to the Monkhorst-Pack scheme [44] is 
generated using the k-point (^, |, ^) with the weight = 1-0 {t-kpoint-rel^ .tme .). The number 
of mesh points of the k-point mesh spanned in the Brillouin-zonc is specified by Lfac. The k-point 
set is then given by the irreducible part of the k-point mesh (c.f. the example in Section TEST 
RUN). 

When the initial wave functions I^E'^k) are calculated by an explicit diagonalization of the Kohn- 
Sham operator in a mixed-basis-set (parameter nbeg) further specifications of the mixed-basis-set 

and of the set up of the initial electron density {init_basis and nrho) are needed. The parameters 
ecuti and tJniLbasis in the file start.inp specify the energy cut-off E"^^ and the atomic pseudo 
orbitals included in the mixed-basis-set. The number of self-consistency cycles of the initialization 
is declared by nomoreJnit. 

The approximation to the exchange-correlation energy functional is determined by the parameter 
Lxc. For applying the GGA correction only in the last step of an LDA calculation set the pa- 
rameter tpostc to . true . . The program applies non-linear core valence exchange and correlation 
automatically, when pseudopotentials that have been generated with a pseudo core density are 
supplied. 

For surface calculations additional contributions to the energy due to a correction for a surface 
dipole moment may be included (parameter tdipol). This requires a surface in the xy-plane [8]. 

The scheme used in the iterative diagonalization is specified by Ledyn. The corresponding parame- 
ters St and 7 of the damped Joannopoulos and William-Soler algorithm are delt and gamma. Note, 
that when the total-energy improvement per iteration is less than eps_chg-delt the parameters delt2 
and gamma2 are used instead, accelerating the convergence of the scheme (with delt>delt2). The 
choice for delt and gamma strongly depends on the atomic species and configuration. However, a 
good guess for 5t lies in the range 1.0 < delt < 20.0 and a good choice for 7 is 7 ~ 0.2. 

In MD simulations and structure optimization the routine f iondyn and f ionsc monitor the error 
of the forces {forcc-eps) before integrating the equations of motion of the nuclei or performing 
a structure optimization step. Note, however, that the forces are not calculated for a specified 
number of iterations {maxjnojorce) and unless the electrons arc sufficiently close to their ground 
state (epse/and epsekinc). The structure optimization terminates, if the residual forces acting on 
the nuclei are sufficiently small (epsfor). The calculation of the electron ground state for a fixed 
configuration of the nuclei is terminated, if the improvement of the total-energy and the wave 
functions per iteration is smaller than epsel and epsekinc rsp. Nevertheless fhi96md is terminated 
on exceeding cither a maximum number of steps (nomore and nstepe) or the overall CPU-time 
limit - which is of importance when the program is running in queuing-system imposing CPU-time 
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restrictions. 



A proper setting of the parameters pfft_store and mesh_accuracy provides an even better perfor- 
mance of the code. The parameter pffLstore determines the fraction of wave functions of which 
the real space representation is stored to avoid a second transformation. The only limitation is 
the available physical memory. The parameter mesh_accuracy specifies the fraction of Fourier co- 
efficients used to represent the electron density (c.f. Appendix B). Choosing mesh_accuracy — 1.0 
implies a proper treatment of the electron density without approximations. In many systems, 
a value of 0.8 results in acceptable loss in accuracy and at the same time in a much better 
performance. 

Output generated during the calculation is written to several files. The chief output file is fort. 6. It 
contains a complete protocol of the initialization, the molecular dynamic simulation or the struc- 
ture optimization and information on each step of the energy minimization. During the molecular 
dynamics simulation the trajectory is also written to the unformatted file fort. 2 (c.f. routine 
f iondyn). The file fort.l contains information on the position of the nuclei and forces, written at 
each self-consistency cycle, while performing structure optimizations or MD simulations. 

Restart files are written every zprmt self-consistency cycles (c.f. inp.mod). The file fort. 71 contains 
among others the wave functions and the coordinates of the nuclei. A restarting run reads this file 
renamed fort. 70. The electron density is stored in file fort. 72. The file fort. 21 contains all necessary 
restart information of the molecular dynamics. 



4 Mciking of the program 

The package fhi96md is available as a tar-archive and can be extracted by the UNIX command tar. 
The directory fhi96md is the root of the package's directory tree. The directory bin contains shell 
scripts for the test run and other examples. These shell-scripts create the input files inp.mod and 
start. inp, compile and run the start utility to generate the input and parameter files for fhi96md and 
finally compile and run the program fhi96nnd. Pseudopotentials for the example runs are included 
in the directory pseudo. They have been generated according to the schemes of Hamann [31]. 
Directories with the generic name example. < scriptname > contain the formatted output of the 
examples and the test run. 

The directory src houses all sources and hbraries of the package. The sources of the utility start 
and the program fhi96md are stored in the corresponding subdirectories. Also included in these 
directories are the makefiles used in the examples to compile the programs by the UNIX command 
make. 

Libraries are contained in the subdirectory lib together with makefiles and sources. These li- 
braries are automatically compiled as recommended by the makefiles of the start utility and the 
program fhi96md. However, there are still some library routines like the fast Fourier transfor- 
mation, ElSPACK-routines and BLAS-routines which we link from commercial libraries such as 
the ESSL. These routines are adopted to the platform running the program and utilize a higher 
performance than public domain routines would do. Currently we offer three versions compatible 
with the ESSL-library, the IMSL-library and the pubhc domain BLAS and FFTPACK routines 
as contained in the netlib. 
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In order to run the program on other platforms, first the makefiles have to be adopted, i.e. the 
FORTRAN?? compiler name as well as the compiler and linker options have to be set properly. 
Automatic zero-initialization of all variables is recommended. 

Second, the available library needs to be adjusted in all example shell scripts as contained in the 
directory bin, i.e. one needs to replace the option essl at the make command call by either imsl 
or netlib. Two routines in libnum have to be ported to the specific platform: the routine cputime 
- measuring the elapsed CPU-time between subsequent calls - and the routine flush - hushing 
the file buffer. 

Note, that some large arrays are declared single precision (real*4 and complex*8) to reduce the 
memory demand. This may require some adjustment in the calls of precision depending routines 
of the external library when single and double precision are used with another convention as e.g. 
on CRAY supercomputers. 



5 Test run 

The test run simply calculates the electron ground state of a bulk GaAs cell. The lattice is a simple 
cubic bravais lattice and the super cell contains eight gallium and arsenic atoms. The calculation 
is performed with an energy cut-off of 8 Ry and 4 special k-points in the irreducible wedge of the 
Brillouin-zone. The initial diagonalization is performed in a basis-set containing plane waves up 
to an energy cut-off of 4 Ry. 

The shell-script run. GaAs. bulk performing the test run creates the input files inp.mod and start. inp, 
compiles and runs the start utility and finally compiles and runs the program fhi96md. The 
workspace is the directory work, where all output of either start and fhi96md is stored. The 
pseudopotential files ga_aa.cpi for gallium and as_aa.cpi for arsenic are provided in the directory 
pseudo and are copied by run. GaAs. bulk to the workspace. The pseudopotentials [40] have been 
created according to the scheme of Hamann [31]. The input files inp.mod and start. inp, the files 
parameter. h and inp.ini, and extracts from the major output file fort. 6 are found at the end of the 
paper. The full set of formatted output files is contained in the directory example. GaAs. bulk. 
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A List of expressions 



In this Appendix we list the central physical quantities and expressions as implemented in the 
routines of the package fhi96md. It is organized in three Subsections: for routines calculating charge 
densities, for routines predominantly dealing with electronic contributions to energy, potentials 
and forces and for routines tabulating structure factors and form factors. The sequence in which 
these routines are called is described in Section 3.1 and in the flow chart in Fig. 3.1. 

Throughout this Appendix and in the package atomic units are used - energies are given in units 
of Hartree - unless noted otherwise and with the exception that in the package reciprocal lattice 
vectors are in units of 27r/aiat- The angular momentum quantum number / translates into an index 
equal to / + 1 in any routine. The weights W]^ of k-points are as given in the input files, though in 
the package the Wk are eventually divided by the number of point group elements nrot. Hence an 
additional factor nrot may appear. The symbols in this Appendix are as defined in the text below 
and the variable name is stated whenever this is of importance. Some general symbols are listed 
in Table A.l. Be aware that a few variables in the package may store different quantities (e.g. in 
rhoe the electron density and the local potential are stored). For the sake of compact formulas, 
summation symbols only bear the summation index, the associated range of summation is listed 
in Table A.2 . 



A.l Charge density of valence electrons and pseudo core 

A. 1.0.1 Routine rhoofr calculates the electron density and the kinetic energy: 
The wave function is given by ^I/i,k(i") = g^^'^ Ui^\s_{r) with: 

;,k(r) = X] ^'^"^ Ci,G+k and j^d^rui,k{r)uj,k{r)^Q5ij, (A.l) 



Ui 



where the ortho-normalization is accomplished by routine graham. 
The electron density n(r) - variable rhoe - is given by 



^(r) = ^EEE^kAkKk(rt)r . (A.2) 

"ski 



The symmetrization is performed as the last step in real space. 
The kinetic energy - variable ekin - is obtained from: 



T' = l E E^kAk E |G + kriQ,G+kP (A.3) 



2 k i G 



A. 1.0. 2 Routine corcha calculates pseudo core charge density and tabulates the form fac- 
tors: 



The pseudocore density n'^°^^{r) is calculated by: 



-core 



» = j:^'^''j:si.{G)^riG) , (a.4) 

G Is 
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where n}^^'^{r) is the form factor of pseudocore density - variable formf-at 

Aqr too 

^r(G) = ^ dTT^U\G\T)nZ'{T) . (A.5) 

The sum in Eq. (A. 4) includes only species for which a pseudo core nj^^^{v) is included in the 
creation of the pseudopotential. 



A. 2 Energy, potential and forces: electronic contribution 

A. 2. 0.3 Routine vof rho calculates the local contributions to energy, potential and forces: 

Due to the long-range tail of the coulomb potential V^{G) and ^/p^^iocaij^Q^ diverge in a periodic 
system, though in the sum of the two potentials the divergent terms cancel. In order to treat the 
potentials separately a charge density n'^^"'^(r) is introduced to remove the divergent contributions 
without affecting the sum of the potentials. The energy contributions E^^'^°^^^ and E™^^~™^^ 
are treated accordingly. By means of n*^^"^(r) these terms are redefined 

For the exphcit expressions see Eqn. (A. 8), (A. 10) and (A. 18). 

The contributions to the local potential in reciprocal space are: 

^locai^Q^ ^ V^{G) + i^p«'i°^=^i(G) + V^^{G) . (A.7) 

In surface calculations an additional contribution \/'^'P°'^(r) to the local potential V^°'^^^{v) arises 
according to Eq. (A. 15) due to a surface dipole moment. 

The Hartree potential is obtained from: 

l^^(G) = -^(n(G)+nG-«(G)) , (A.8) 



n = 1 / dVe-'^"n(r) (A.9) 



with 

n(G) 
and 

^Gaufi(Q) = J2 SiXG)^f;''\G) , (A.IO) 
Is 

where $f'^'^^(G) is given by Eq. (A. 35). 
The local pseudopotential V^'^'^°^^^{G) reads: 

yps,local(Q) _ j2 Si^{G)<^f^{G) (A.ll) 
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with <S'/^(G) and $j^^(G) given by Eqn. (A. 33) and (A. 34) respectively. 
The exchange-correlation potential is calculated from: 

where the pscudocorc density is inchidcd, if provided with the pscudopotcntials. is ap- 

proximated by LDA, BP, or PW91 functional and in the latter two cases the first and second 
derivatives of the electron and pseudocore density are evaluated in routine sivxc^gc: 



Vn(r) = ^iGe^^-'"n(G) (A.13) 

G 

92 



-n\ 



r) = -EG'Ae'G"^(G) ■ (A.14) 



The potential due to a surface dipole in an orthorhombic cell with the surface in the a;|/-plane, 
i.e. lattice vector as = 0362, is given by: 

^dipole^r) = [z e{z^^ -Z) + {Z- Z"^^) 9{Z - Z^'^))E'''''' + V^,,dipole , (A.15) 

with the surface dipole moment d according to Eq. (A. 16) and where z™"^ is the coordinate at 
which n^'^^{z) — J dxdyn{v) has its minimum. The surface dipole moment is given by 

d = d'^- , (A. 16) 

where the electric dipole is 



i-z'—'-i-aa ,■ 

d^ — dz dxdyn{r) 

y^min J 



and the ionic dipole is 

c?''"^ = E ('^^s./a • e. - z'^''' + 9{z^'^ - Ti^j^ ■ e,) - as 



with H.dipole 



and where E"^^'^ is 



fi Id f ■ 
H.dipole — —E ^ ( (03 — ^™°) 



£;fieid^_l^^ _ (A.17) 



The Energy expressions evaluated in vof rho read: 

^kin _|_ _|_ ^sr ^self _j_ ^ps, local 

+£;PS,nl ^ ^XC ^ ^dipole {k.lS) 
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and 



^Harris ^ ^ ^ ^ ^.(k) _ E'^'^'Iu^^-'^] + E^'V-'^] 

k i 

- V^^[n^^-^^] + EH,Gaufi ^ ^sr _ ^self _ ^dipole ^^_^g^) 

with n^'^~^^(r) being the charge density of the previous iteration and T^, E^^, E^^^^ and E^^ after 
Eqn. (A.3), (A.39), (A.36) and (A.25) rsp. 
The Hartree energy is given by: 

G^O 1^1^ 

The local pseudopotential energy is obtained from: 



^ps,local ^ Q ^ yps,locaI(Q) ^(q) (^ 21) 

G 

The exchange-correlation energy reads: 

£;xc = f (n{r) + n"°"'^(r)) (n(r) + ^""'^^(r)) (A.22) 
Jn 

and the exchange-correlation potential energy is 

yxc ^ ^ dhn{r) (n(r) + n^°'^^(r)) , (A.23) 

where e^*^ (n(r)) is approximated by LDA, BP, or PW91 (c.f. V'^^ Eq. (A.12)). 
Energy of the surface dipole E^'^'P"'*^ is: 

^dipol ^ _£'field ^ 

with the dipole moment d and the electrostatic field E^^^^ according to Eqn. (A. 16) and (A. 17) 
rsp. 

The contributions to Harris energy ^jHams ^ead: 



G^O 

^H,Gaufi ^ 27r O ^ ■ 
G^O 



|n(G)p 
|G|2 

^Gaufi/QN|2 



IG 



2 



The local and ionic contribution to forces - variable fion - are given by: 
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F/.,/. = - iO ^ Ge'^- (n(G) + nGaufi(G))$J-«(G) 



G 



+$P:(G)n(G) (A.24) 



with according to Eq. (A. 41). 



A. 2. 0.4 Routine nlrhkb calculates the non-local contributions to energy, potential and forces: 
The non-local contribution to energy E^^ - variable enl - is given by: 

^"'=EEEE^kAk<j/^:L,.,,^(k)i^ , (A.25) 

Inja k l,m i 

with w^^i and flj^j^^i^rni^) ^-fter Eq. (A. 38) and Eq. (A. 28) rsp. If the nuclei are located on a 
mesh commensurate with the super cell the numerical effort in evaluating (A.25) can be reduced 
significantly as described in [35] - see also Appendix B; the necessary phase factors and other 
relevant quantities are tabulated by routine phf ac. 

The non-local contribution to forces Fj^j^ -variable fiori-nl- are computed from: 

Ki. = EEE^kAk<,Im{f,,,,,,,,^(k)i-)^,,^,,^(k)} , (A.26) 

where f^Xj^^i^rnO^) and fi,i,,i,,i,mO<-) -variables fnl, sfl, sf2 and sfS- are given by: 

f^,i.,i,M = E e-^^^^-^^--^'^^ $?:;i(G + k) c,G+k (A.27) 

G 

and 

UjMmO^) = If E(G + k)e-'(^+^)-=.^^$?::i(G + k) c,G+k , (A.28) 
with ^f^'^rni^ + k) according to Eq. (A.37). 

A. 2. 0.5 Routine df orce computes the application of Hamiltonian to wave function: 

(G + k\H^^\^i^k) = (G + k|f + t>^°"^' + t>P''°^|^i,k) . (A.29) 
The kinetic contribution reads: 

(G + k|f|^i,k) = ^|G + k|2c,,G+k . (A.30) 
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The local potential contribution (obtained by FFT) is given by: 

(G + k|y'°-i|*,,k) = ^ dr' y'°^^(r) Ui,^{r) e-^^"" . (A.31) 
The non-local potential contribution is obtained from: 

(G + k|VP^'-'|^,0 = E E KJSmM e'^-^^-^^^?::i(G + k) (A.32) 

IsJa, l,m 

with w'^li and $/')°^(G + k) according to (A.38) and (A.37) rsp. 

The same optimization as in routine nlrhkb for nuclei located on a mesh commensurate with the 
super cell applies here as well. 

A. 3 Structure and form factors and ionic contributions to energy and forces 

A. 3. 0.6 Routine struct tabulates the structure factor of the ionic basis -variable struct: 

5,,(G) = Ee^^-^^-^^ . (A.33) 



A. 3. 0.7 Routine f ormf tabulates the form factors of Gaussian charges and of the pseudopo- 
tential and calculates the electrostatic self-energy of the Gaussian charges: 

The form factor of the local potential -variable vps- is given by: 

^ r + T^erf } . (A.34) 
The form factor of GauB charges -variable rhops- reads: 

The electrostatic self-energy of Gaussian charges -variable eself- is calculated from: 

E-'' = 4tr,nfr . (A.36) 

V is is 



A. 3. 0.8 Routine nlskb tabulates the form factor and the prefactor of non-local pseudopo- 
tentials: 

The form factor of non-local pseudopotentials -variable pkg, pkg_a- are obtained from: 
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where 



and yim{dG,(f>G) are 



,m> 

,171 — 

, m < 



(A.37) 



The prefactor wj^i -variable wnl- reads: 



ft 



-1 



(A.38) 



A. 3. 0.9 Routine ewald calculates the screened ionic contributions to the energy and the 
forces: 



The energy of screened ions -variable esr- is given by: 



R| 



erfc 



,Gaufi2 , ^Gaufi2 



(A.39) 



where in the innermost sum (0, Jg, Jg) ^ (R, Is, la) is obeyed. The forces due to screened nuclei- 
nuclei interaction -variable fionsr- is obtained from: 



= E E E <ii.<ij. 



\ri,j, - Tj^,j^ - Rp 



1 Irr / — T/ 7 — R| 

X < — I exp 

- rj^,j, - R| 



|tJ3,J^ - Tj^,J^ - R| 
„Gaufi2 I Gaufi2 



+ erfc 



,Gaufi2 I Gaufi2 
r. "T ' J. 



(A.40) 



where in the innermost sum (0, Jg, Jg) 7^ (R, 4; -^a) is obeyed. The sum over R is truncated at 
large R''''\ 



B The parameter file parameter. h and the input file inp.ini 



The parameter file parameter. h and the input file inp.ini are usually generated by the start utility 
start from the files inp.nnod and start.ini. The program fhi96md runs also individually without the 
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help of the start utihty start. This requires the user to provide the files parameter. h and inp.ini in 
addition to the pseudopotential files and the optional restart files. 



In Tables B.l and B.2 we describe the parameter file parameter. h and the input file inp.ini. It 
should be noted that the file inp.mod contains the entries mesh-accuracy and pfft_store, which 
have no effect in the program fhi96md and are required only by the utility start to generate the 
parameters nfftstore and nrlx, nr2x, and nr3x. 

The file para meter, h is included by almost all source files at compilation time and is a requisite 
to compile the program fhi96md. Consequently this file has to meet FORTRAN?? syntax rules 
and deviations from these may result in errors at compilation time. All of the parameters listed 
in Tab. B.l have to be declared as integer variable and FORTRAN?? parameter statements are 
used to assign the corresponding values as described in this Table. The parameter n^'wx determines 
the maximum number of plane waves used to represent the wave function. For a given energy 
cut-off ii^cut (in Ry) it should be set according to 

1 3 

ngwx > ^ ^ -E'cut ■ 



ngwix for a given i?™** should be set accordingly. The size of the Fourier mesh determines the 
accuracy of the charge density. Set the parameters nrlx, nr2x and nrSx according to the sampling 
theorem: 



nrlx > - ||ai|| V^^cut 



and nr2x, nrSx correspondingly. Using smaller values for nrlx, nr2x and nrSx means to skip the 
highest G-vectors in 

^(r) = EEEE^.kC-^^.G-,keH«-«)- 

k i G G 



and results in a better performance. However, the applicability of the grid should be checked for 
each system individually. In particular in systems with strongly localized orbitals this may be an 
unacceptable approximation. 

The entries in file inp.ini are hsted in Tab. B.2. When separated by comma, they are expected on 
the same line of file inp.ini. Whenever the entries have the same meaning as in the file start. inp, 
we refer the reader to Tab. 3.2. 

For each species na lines containing coordinates and the boolean parameters tford and t^auto-coor 
are expected. Lines containing the velocity of a nucleus should immediately follow the line with 
the coordinates of the corresponding ion. Whenever some ions of a species are located on a mesh 
commensurate with the super cell, the evaluation of the non-local contributions to energy and 
potential may be accomplished in a more efficient fashion. This requires the variable tford of the 
relevant ions to be set to .false, and the variable t^auto-coord to be set to .true.. The entry ineq^pos 
contains the number of mesh points per super cell along each lattice vector. Note, that the point 
(0,0,0) in the super cell must be a point on the mesh. 



25 



Acknowledgements 

The authors are indebted to Dr. E. Pehlke for many valuable discussions during the development 
of this version of the package. 

References 

[I] W. Andreoni, F. Gygi, and M. Parrinello, Phys. Rev. Lett. 68, 823 (1992). 
[2] N. Troullier and J. L. Martins, Phys. Rev. B 46, 1754 (1992). 

[3] G. Ortiz, Phys. Rev. B 45, 11328 (1992). 

[4] A. Garcia et al, Phys. Rev. B 46, 9829 (1992). 

[5] J. D^browski and M. Scheffler, Phys. Rev. B 40, 10391 (1989). 

[6] S. B. Zhang and J. E. Northrup, Phys. Rev. Lett. 67, 2339 (1991). 

[7] J. A. Alves, J. Hebcnstrcit, and M. Scheffler, Phys. Rev. B 44, 6188 (1991). 

[8] J. Neugebauer and M. Scheffler, Phys. Rev. B 46, 16067 (1992). 

[9] R. Stumpf and M. Scheffler, Phys. Rev. B 53, 4958 (1996). 

[10] M. J. Caldas, J. D^browski, and M. Scheffler, Phys. Rev. Lett. 65, 2046 (1990). 

[II] E. Pehlke and M. Scheffler, Phys. Rev. Lett. 71, 2338 (1993). 
[12] O. Pankratov and M. Scheffler, Phys. Rev. Lett. 75, 701 (1995). 
[13] M. Bockstedte and M. Scheffler, Z. Phys. Chem. in print. 

[14] A. Grofi, B. Hammer, M. Scheffler, and W. Brenig, Phys. Rev. Lett. 73, 3121 (1994). 
[15] E. Pehlke and M. Scheffler, Phys. Rev. Lett. 74, 952 (1995). 
[16] C. Stampfl and M. Scheffler, Phys. Rev. Lett. 78, 1500, (1996). 
[17] N. Moh et al, Phys. Rev. B 52, 2550 (1995). 

[18] A. Grofi, M. Bockstedte, and M. Scheffler, submitted to Phys. Rev. Lett. 

[19] O. Sugino and R. Car, Phys. Rev. Lett. 74, 1823 (1995). 

[20] P. Hohenbcrg and W. Kohn, Phys. Rev. 136B, 864 (1964). 

[21] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 567 (1980). 

[22] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 

[23] A. D. Becke, Phys. Rev. A 38, 3098 (1988). 

[24] J. P. Perdew, Phys. Rev. B 33, 8822 (1986). 

[25] J. P. Perdew et al, Phys. Rev. B 46, 6671 (1992). 

[26] W. Kohn and J. L. Sham, Phys. Rev. 140A, 1133 (1965). 



26 



[27] J. Ihm, A. Zunger, and M. L. Cohen, J. Phys. C 12, 4409 (1979). 
[28] A. Williams and J. Soler, Bull. Am. Phys. Soc. 32, 562 (1987). 
[29] M. C. Payne et al, Phys. Rev. Lett. 56, 2656 (1986). 

[30] G. B. Bachelet, D. R. Hamann, and M. Schliiter, Phys. Rev. B 26, 4199 (1982). 

[31] D. R. Hamann, Phys. Rev. B 40, 2980 (1989). 

[32] N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991). 

[33] X. Gonze, R. Stumpf, and M. Scheffler, Phys. Rev. B 44, 8503 (1991). 

[34] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982). 

[35] R. Stumpf and M. Scheffler, Comp. Phys. Comm. 79, 447 (1994). 

[36] S. Nose, J. Chem. Phys. 81, 511 (1984). 

[37] W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 

[38] J. Hebenstreit, M. Heinemann, and M. Scheffler, Phys. Rev. Lett. 67, 1031 (1991). 
[39] S. G. Louie, S. Proyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982). 
[40] M. Fuchs and M. Scheffler, in preparation. 

[41] M. Weinert and J. W. Davenport, Phys. Rev. B 45, 13709 (1992). 

[42] M. J. Gillan, J. Phys. Cond. Matt. 1, 689 (1989). 

[43] F. Wagner, T. Laloyaux, and M. Scheffler, submitted to Phys. Rev. B. 

[44] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 

[45] D. J. Chadi and M. Cohen, Phys. Rev. B 8, 5747 (1973). 

[46] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

[47] F. Tassone, F. Mauri, and R. Car, Phys. Rev. B 50, 10561 (1994). 

[48] M. C. Payne et al, Rev. Mod. Phys. 64, 1045 (1992). 

[49] J. Ncugcbauer and C. G. van dc Walle, in Materials Theory, Simulations, and Parallel Algorithms, 
Vol. 408 of MRS Symposia Proceedings, edited by E. Kaxiras, J. Joannopoulos, P. Vashita, and R. K. 
Kaha (MRS, Pittsburgh, Pennsylvannia, 1995). 

[50] D. Singh, H. Krakaucr, and C. Wang, Phys. Rev. B 34, 8391 (1986). 

[51] L. Verlet, Phys. Rev. 159, 98 (1967). 

[52] G. Pastore, E. Smargiassi, and F. Buda, Phys. Rev. A 44, 6334 (1991). 
[53] P. E. Blochl and M. Parrinello, Phys. Rev. B 45, 9413 (1992). 

[54] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Claredon Press, Oxford, 1990). 

[55] J. Stoer and R. Bulirsch, Einfiihrung in die Numerische Mathematik, 2 ed. (Springer- Verlag, Berlin, 
1978), Vol. II. 

[56] M. R. Pederson and K. A. Jackson, Phys. Rev. B 43, 7312 (1991). 



27 



TEST RUN 



Include file inp.mod 



-1 100 1000000 
100 1 
12.0 0.2 

4.0 0.2 0.0001 
400 2 
0.0 1.0 
2 2 

.false. 

.F. 0.001 .F. 0.002 

.false, .false, .false. 1800 

.false. 

0.0001 0.0005 0.1 

0.01 0.01 3 

1 



nbeg iprint timequeue 
nomore nomore_init 
delt gamma 

delt2 gamma2 eps_chg_dlt 
delt_ion nOrder 
pfft_store mesh_accuracy 
idyn i_edyn 
i_xc t_postc 
trane ampre tranp amprp 
tfor tdyn tsdp nstepe 
tdipol 

epsel epsfor epsekinc 
force_eps mcLX_no_f orce 
init_basis 



Include file start. inp 



2 

5 

1 

10.47 0.0 0.0 
1 

0.5 0.5 0.5 1.0 

3 3 3 
.false . 
8 4.0 

0.004 .true. .f. 

.true, .false. 1 

5 2 1234 

873.0 1400.0 le8 1 

.t. .true. 
4 3 'Gallium' 
1.0 3.0 0.7 3 3 
.t. .t. .f. 
0.0 0.0 0.0 .t. 
0.5 0.5 0.0 .t. 
0.5 0.0 0.5 .t. 
0.0 0.5 0.5 .t. 
4 5 'Arsenic' 
1.0 3.0 0.7 3 3 



number of species (nsp) 
excess electrons 
number of empty states 
ibrav pgind 
celldm 

number of k-points 

k-point coordintes, weight 

fold parameter 

t_kpoint_rel 

Ecut [Ry] , Ecuti [Ry] 

ekt tmetal tdegen 

tmold tband nrho 

npos nthm nseed 

T_ion T_init Q nfi_rescale 

tpsmesh coordwave 

number of atoms, zv, name 

gauss radius, mass, damping, l_max, l_loc 

t_init_basis 

tauO tford 

tauO tford 

tauO tford 

tauO tford 

number of atoms, zv, name 

gauss radius, mass, damping, l_max,l_loc 
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.t. .t. .f. 

0.25 0.25 0.25 
0.75 0.25 0.75 
0.75 0.75 0.25 
0.25 0.75 0.75 



: t_init_basis 

.t. : tauO tford 

.t. : tauO tford 

.t. : tauO tford 

.t. : tauO tford 



Include file para meter, h - generated by the start utility 



c========= Parameters for cp : ================ 

integer nsx , nax , nx , ngwx , ngx 

integer ngwix , nx_init , nr Ix , nr2x , nrSx , nschltz 

integer nx_basis ,max_basis_n,nlmax_init 

integer nnrx , nkptx , nlmsLX , mmaxx , n_f f t_st ore 
c SPECIES ATOMS STATES 

parameter (nsx= 2, nax= 4, nx= 21) 
c FULL BASIS 

parameter (ngwx= 449, ngx= 3592) 
c INITIAL BASIS 

parameter (ngwix= 168, nx_init= 169) 

parameter (nx_basis= 1, max_basis_n= 21) 

parameter (nlmax_init= 1) 
c FFT: X-MESH Y-MESH Z-MESH 

parameter (nrlx= 20,nr2x= 20,nr3x= 20) 

parameter (nnrx= 8400) 
c K-POINTS 

parameter (nkptx= 4) 
c Number of ffts to be stored between rhoofr and dforce 

parameter (n_fft_store= 1) 
c Number of Im-components and max length of radial mesh 

parameter (nlmax= 4,mmaxx= 550) 

parameter ( nschltz = 1 ) 
0========= end of parameters for cp : ========== 



Include file inp.ini - generated by the start utility 



1 :ibrav, pgind 

32.0000 T .00400 F : nel,tmetal,ekt,tdegen 

8.00000 4.00000 : ecut,ecuti 

T F 1 : tmold,tband,nrho 

5 2 1234 : npos, nthm, nseed 

873.00 1400.00 . lOOOE+09 1 : T_ion, T_init, Q, nfi.rescale 

2 T T :nsp,tpsmesh,coordwave 
4 : nkpt 

.1666667 .1666667 .1666667 .2962963 :xk(l-3) ,wkpt 
.1666667 .1666667 .5000000 .4444444 :xk(l-3) ,wkpt 



29 



. 1666667 
.5000000 
10.47000000 
.00000000 
.00000000 
1.00000000 
.00000000 
.00000000 
10.4700000 
'Gallium ' ^ 
.70000 1.00000 
T T F 
.000000000 
5.235000000 
5.235000000 
.000000000 

'Arsenic ' ^ 
.70000 1.00000 
T T F 
2.617500000 
7.852500000 
7.852500000 
2.617500000 



24 
1 



5000000 
5000000 

.00000000 
10.47000000 
.00000000 
. 00000000 
1 . 00000000 
.00000000 

1147.73082300 
3.00000 3.00000 



5000000 . 2222222 
5000000 . 0370370 
.00000000 
.00000000 
10.47000000 
.00000000 
.00000000 
1.00000000 



xk(l-3) ,wkpt 
xk(l-3) ,wkpt 
lattice vector al 
lattice vector a2 
lattice vector a3 
rec. lattice vector bl 
rec. lattice vector b2 
rec. lattice vector b3 
: al at, omega 
name, number, valence charge, ion_fac 



3 3 



: ion_dainp,rgauss,l_max,l_loc 
s,p,d 

.000000000 

.000000000 
5.235000000 
5.235000000 



: t_init_basis 
.000000000 
5.235000000 
.000000000 
5.235000000 
; ineq_pos 

5.00000 3.00000 : name, number , valence charge, ion_fac 
3 3 : ion_damp,rgauss,l_max,l_loc 

:t_init_basis s,p,d 



F 

F 
F 
F 



F 

F 
F 
F 



T 
T 

T 
T 



2.617500000 
2.617500000 
7.852500000 
7.852500000 



2.617500000 F F F T 

7.852500000 F F F T 

2.617500000 F F F T 

7.852500000 F F F T 
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: ineq_pos 








nrot = number 


1 













1 













1 




1 
















1 







1 







1 













-1 













1 





Output file file fort. 6 

******* this is the complex fhi96md program ******** 

******* ibm - version ******** 

******* Juli 1996 ******** 

»>nbeg= -1 nomore= 100 iprint= 100 
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> Exchange: LDA 

>============================================ 

»>electronic time step= 12.0000 gaiiima= .2000 

»> Using delt= 4.0000 gaiiima= .2000 when energy changes less than: 
>accuracy for convergency: epsel= .00010 epsfor= .00050 epsekinc= .10000 
>damped Joannopoulos algorithm for electron dynamics 
>ions are not allowed to move 
normally no mixing of old charge is done 
>ibrav= 1 pgind= nrot= 24 alat= 10.470 omega= 1147.7308 mesh= 20 20 20 
>ecut= 8.0 ryd ecuti= 4.0 ryd nkpt= 4 

nel, tmetal, ekt, tdegen= 32.0000000000000000 T 0.400000000000000008E-02 F 
>alat= 10.470000 alat= 10.470000 omega= 1147.730823 

lattice vectors 
al 10.470000 .000000 
a2 .000000 10.470000 

a3 .000000 .000000 

reciprocal lattice vectors 



lOOOE-03 



bl 

b2 
b3 



1 . 000000 

.000000 
.000000 



. 000000 
1.000000 
.000000 
positions tauO 
specie Nr. x y 

Gallium 1 .0000 .0000 

Gallium 2 5.2350 5.2350 
Gallium 3 5.2350 .0000 5 

Gallium 4 .0000 5.2350 5 

Arsenic 1 2.6175 2.6175 2 
Arsenic 2 7.8525 2.6175 7 
Arsenic 3 7.8525 7.8525 2 
Arsenic 4 2.6175 7.8525 7 
nkpt= 4 

weight of all kpts: 0.999999900000000053 
... so I'll scale them for you... 
ratios of FFT mesh dimensions to sampling theorem 
>ps-pots as given on radial mesh are used 
>gvk: ngwx and max nr. of plane waves: 449 



.000000 
.000000 
10.470000 

.000000 

.000000 
1.000000 



z 

,0000 

,0000 
,2350 
,2350 
,6175 
,8525 
,6175 
,8525 



1.061 1.061 1.061 



440 



k-point weight # of g-vectors 

1 .17 .17 .17 .0123 440 

2 .17 .17 .50 .0185 434 

3 .17 .50 .50 .0093 432 

4 .50 .50 .50 .0015 432 
Weigthed number of plane waves npw: 435.248 
Ratio of actual nr. of PWs to ideal nr.: .99246 

# of electrons= 32.0000, # of valencestates= 16, # of conduction states= 5 

atomic data for 2 atomic species 

pseudopotentialparameters for Gallium 
>nr. of atoms: 4, valence charge: 3. 000, force fac: 
l_max 3 l_loc: 3 rad. of gaussian charge: 1.000 



3.00, speed damp: .7000 
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pseudopotentialparameters for Arsenic 

>nr. of atoms: 4, valence charge: 5. 000, force fac: 3.00, speed damp: .7000 
l_max 3 l_loc: 3 rad. of gaussian charge: 1.000 
Final starting positions: 

.0000 .0000 .0000 5.2350 5.2350 .0000 5.2350 .0000 5.2350 

.0000 5.2350 5.2350 2.6175 2.6175 2.6175 7.8525 2.6175 7.8525 

7.8525 7.8525 2.6175 2.6175 7.8525 7.8525 



phfac: is, n_ideal: 1 
phfac: is, n_ideal: 2 

phfac : is , i_kgv,n_class 10 

phfac: is, i_kgv,n_class 2 
>nlskb:is= 1 wnl: 1 -.1067262 2 -.7541519 3 -.7541519 4 -.7541519 
>nlskb:is= 2 wnl: 1 -.1008436 2 -.4304406 3 -.4304406 4 -.4304406 



starting density calculated from pseudo-atom 



formfa: rho of atom in 3.000000 
formfa: rho of atom in 5.000000 



SYMMETRY OPERATIONS 



>s(isym) in latt . 




>sym( . . ) in cart . 

.00000 .00000 

>s(isym) in latt. 
1 



coord: 1 
1 

coord: 1.00000 
.00000 1.00000 
coord: 1 
1 





.00000 .00000 




1 

. 00000 1 . 00000 




CHECK SYMMETRIES 



>Center of symmetry symO= .000000 .000000 

Table of symmetry relations of atoms, 
symmetric at., xneu=tauO+<sym.Op. (isym)>) 

is ia iasym isym tauO xneu 

1 1 1 1 .00000 .00000 .00000 
. 00000 

11 12 
. 00000 

2 4 2 23 
7.85250 

2 4 2 24 
7.85250 



.00000 .00000 



.00000 



2.61750 7.85250 7.85250 
2.61750 7.85250 7.85250 



.000000 

(iasjrm=nr. of 



,00000 .00000 



.00000 .00000 



-2.61750 -7.85250 



-2.61750 -7.85250 
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ITERATIONS IN INIT STARTED 



ITERATIONS IN MAIN STARTED 



=== LOOP n_it= 1 
time elapsed for nlrh t = .2900 
rhoofr stores starting density for mixing in c_fft_store 

rhoofr: integrated electronic density in g-space = 31.999267 in r-space = 
31.999255 

time elapsed for rhoofr t = 1.0000 



internal energy at zero temperature 




-33, 


.989843 


a. 


,u 


non-equillibrium ( 


antropy 






.000000 


kB 


equillibrium ( 


sntropy 






.000000 


kB 


kT 


energy 






.004 eV 






(non-eq) free 


energy 




-33, 


. 989843 


a. 


.u 


(non-eq) total 


energy 




-33, 


. 989843 


a. 


,u 


Harris 


energy 




-34, 


.290291 


a. 


,u 


kinetic 


energy 




11, 


.557750 


a. 


,u 


electrostatic 


energy 




-40. 


.196264 


a. 


,u 


real hartree 


energy 




2, 


.813098 


a. 


,u 


pseudopotential 


energy 




6, 


.191265 


a. 


,u 


n-1 pseudopotential 


energy 




-1, 


.989500 


a. 


,u 


exchange-correlation 


energy 




-9. 


.553094 


a. 


.u 


exchange-correlation potential 


energy 




-12, 


.456593 


a. 


,u 


kohn - sham orbital 


energy 




-2, 


.810464 


a. 


, u 


self 


energy 




54, 


.256150 


a. 


,u 


esr 


energy 






. 000307 


a. 


,u 


gaussian 


energy 




22, 


.685616 


a. 


.u 



time elapsed for n x nkpt x graham/ortho = .1200 

nel dampeig true_efermi efermi ekt seq sneq 
32.0000 .700 3.20000 3.20000 .004 .00000 .00000 

k-point .167 .167 .167, eigenvalues and occupation numbers: 
eig -9.053 -7.144 -7.144 -7.144 -3.100 -3.100 -3.100 -.402 -.402 -.402 
eig .557 .577 .577 .577 2.885 2.885 4.859 4.859 4.859 5.076 
eig 6.576 
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occ 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
occ 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 .0000 .0000 .0000 .0000 
occ .0000 



>»n_it nfi Ekinc Etot Eharr Ezero mPorce mChange 

Seq Sneq Efermi Dvolt W_a 
»> 1 4.67996 -33.98984 -34.29029 -33.98984 .00000 .000 
.0000 .0000 3.2000 .0000 .0000 

»> OK, I stop after timestep nr. 10000 
time elapsed per electronic time step t = 3.3800 
time in queue: 1000000 max. number of steps: 281060 

»>fermi: No. of eigv with eig_force > 10% ekt: 84 
»> 2 7.31348 -34.26236 -33.98984 -34.26236 .00000 .000 
.0000 .0000 3.2000 .0000 .0000 



»> 18 .06177 -34.35555 -34.35555 -34.35555 .00000 .000 
.0000 .0000 3.2000 .0000 .0000 
=== LOOP n_it= 19 
phfac: is, n_ideal: 2 

rhoofr: integrated electronic density in g-space = 32.000000 in r-space = 
31.999986 



internal energy at zero temperature 




-34, 


.355551 


a. 


,u 


non-equillibrium ( 


antropy 






.000000 


kB 


equillibrium ( 


sntropy 






,000000 


kB 


kT 


energy 






,004 eV 






(non-eq) free 


energy 




-34, 


,355551 


a. 


,u 


(non-eq) total 


energy 




-34. 


,355551 


a. 


,u 


Harris 


energy 




-34. 


.355552 


a. 


, u 


kinetic 


energy 




11. 


,920575 


a. 


,u 


electrostatic 


energy 




-40, 


. 142623 


a. 


,u 


real hartree 


energy 




2, 


.909120 


a. 


,u 


pseudopotential 


energy 




5. 


.960148 


a. 


,u 


n-1 pseudopotential 


energy 




-2. 


.496765 


a. 


,u 


exchange-correlation 


energy 




-9. 


.596886 


a. 


• U 


exchange-correlation potential 


energy 




-12. 


.514237 


a. 


.u 


kohn - sham orbital 


energy 




-2. 


.793552 


a. 


.u 


self 


energy 




54. 


.256150 


a. 


,u 


esr 


energy 






.000307 


a. 


,u 


gaussian 


energy 




22. 


.685616 


a. 


,u 



&&s atomic positions and local+nl forces on ions: 
Gallium : 

>&&s-n .000000 .000000 .000000 .000000 .000000 

.000000 

>&&s-n 5.235000 5.235000 .000000 .000000 .000000 .000000 
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>&&s-n 
>&&s-n 



5 . 235000 
.000000 



.000000 
5.235000 



5 . 235000 
5.235000 



,000000 
,000000 



.000000 
.000000 



.000000 
.000000 



Arsenic : 

>&&s-n 2.617500 2.617500 2.617500 

>&&s-n 7.852500 2.617500 7.852500 

>&&s-n 7.852500 7.852500 2.617500 

>&&s-n 2.617500 7.852500 7.852500 
>simi of all (local+nl) forces / n_atoms = 

. 0000000000 
(should = 0) 

nel dampeig true_efermi efermi 

32.0000 .700 3.20000 3.20000 



. 000000 
. 000000 
.000000 
.000000 
,0000000000 



.000000 
.000000 
.000000 
.000000 

.0000000000 



.000000 
.000000 

.000000 
.000000 



ekt seq sneq 
.004 .00000 .00000 



> 1. k-point .167 .167 .167, ngw 440, EWs and OCCs: 
>eig: -9.370 -7.429 -7.429 -7.429 -3.427 -3.427 -3.427 -.549 -.549 
-.549 .260 .438 .438 .438 2.537 

>eig: 2.537 4.705 4.705 4.705 4.847 6.349 

>occ: 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
>occ: 2.0000 .0000 .0000 .0000 .0000 .0000 

> 2. k-point .167 .167 .500, ngw 434, 



-8.621 
.318 
1.277 
2 . 0000 



>eig: 
-.955 
>eig: 
>occ : 
2.0000 2.0000 
>occ: 2.0000 
k-point 



> 3. 
>eig: 
-.039 
>eig: 
>occ : 



-8.621 -7.350 -7.350 -3.535 
.318 .997 .997 1.277 
5.550 5.550 6.761 6.761 6.957 
2.0000 2.0000 2.0000 2.0000 2.0000 2 
2.0000 2.0000 2.0000 2.0000 
.0000 .0000 .0000 .0000 .0000 
.167 .500 .500, ngw 432, 



EWs and OCCs: 
3.535 -1.754 -1.754 -.955 



0000 2.0000 2.0000 



EWs and OCCs: 



-8.033 -8.032 -8.032 -8.032 -3.061 -3.061 -3.061 -3.061 -.039 
-.039 -.039 1.684 1.684 1.684 
1.684 5.608 5.608 5.608 5.608 7.767 
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
>occ: 2.0000 .0000 .0000 .0000 .0000 .0000 

> 4. k-point .500 .500 .500, ngw 432, EWs and OCCs: 

>eig: -7.951 -7.951 -7.951 -7.951 -3.815 -3.815 -3.815 -3.815 1.756 
1.756 1.756 1.756 1.756 1.756 1.756 
>eig: 1.756 4.192 4.192 4.192 4.192 7.634 

>occ: 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 
2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 



>occ: 2.0000 .0000 .0000 
>»n_it nfi Ekinc Etot 
Seq Sneq Efermi 



0000 .0000 
Eharr 

Dvolt W_a 



»> 



19 
,0000 



av. 



.05767 -34.35555 -34.35555 
.0000 3.2000 .0000 .0000 
=========== END OF THE MAIN-LOOP ======= 

time elapsed for nlrh t = .2911 



.0000 
Ezero 

-34.35555 



mPorce 



mChange 



.00000 



.000 
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av. time elapsed for rhoofr t = .9821 

av. time elapsed for vofrho t = .2632 

av. time elapsed for n x nkpt x dforce = 1.4511 

av. time elapsed for nkpt x graham/ortho = .1063 

av. time elapsed for rest (in main) t = .1979 

av. time elapsed per elec. time step t = 3.2916 
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Figure captions: Fig 3.1:Flowchart of the program fhi96md. The routines init and f iondyn 

are described in more detail in the text. Output is generated at 
the end of each self-consistency cycle and by the routines f iondyn, 
f ionsc, fermi, init and vofrho. Restart files are written by the 
routine f iondyn and by a call of routine o_wave in the main pro- 
gram. 
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initializat 


on: 


header: 


read inp.mod 


init: 


read inp.ini and pseudo potentials 




calculate structure and form factors 




(f ormf , nlskb, phf ac, struct) 




set up initial electron density and 




wave functions 



calculate 




rhoof r: 


electron density n^'^\r) 


nlrhkb: 


non-local contributions 




to energy, potential (and forces) 


vof rho: 


local contributions 




to energy, potential (and forces) 




move atoms (structure optimization or dynamics) 



compute 

dforce: (G + k\nKs[n^'\r)]\¥:\\i}) 
iterate wave functions |^'|j^^(k)) 
l*M(k)) ^ |f l^+^'(k)) 



repeat for all states 



perform 






graham: Gralim 


-Schmidt ortlionormalisation 




k)V.. 






repeat for all k-points 


^ 



r ^ T + 1 



calculate (e.g. for metals): 

f ermi: occupation according to fermi 
distribution 

validate 

iteration parameters (c.f. files delt, 
ion_fac, stopfile and stopprogram). 
I 

'terminate 

on convergence 
or 

if number of iterations or 
cpu-time is exceeded 

terminate 



call 

f ionsc: for structure optimization 

f iondyn: molecular dynamics 
both routines monitor the error of the forces and 
move the ions after convergence. 

\~ 

have atoms been moved ? / 



yes 



recalculate 




struct. 




phf ac: 


structure and phase factors 


calculate 




nlrhkb: 


non-local contributions 

to energy, potential (and forces) 


vof rho: 


ewald sum and local contributions 
to energy, potential (and forces) 



Fig. 3.1. 
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Table 3.1. input file inp.mod 
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inp.mod 



T-\ Q T» Q m fi'i" dT 


ylV^I rdiige 




nbeg 




set up of the initial wavefunction: 




-1 


■0init by an diagonalization on a subset of the plane wave 
uasis-ser (^c.i. ioec. z.z j 




-2 


^init is read from the file fort. 70 


iprint 


integer 


number of electronic iterations between detailed output and 
writing oi resxarx nies yz.L. texx j 


timequeue 


integer 


maximum CPU-time in seconds: 

the program terminates before exceeding this time and 
restart files are written 


nomore 


integer 


maximum number of : 

11 ^^yor, zo/yn — . raise. j. seii-consistency cycies 
11 \oj07 or Ziiyii — .Xiucj. cLTjOiiiic niovcb 


liOlUOl 6-lllli 


iiitegei 


iiiaxiiiiuiii iiLiiiiijei Oi oeii-coiibioieiicy cycieo iii uiie iiiiiiaiioa- 
tion 


delt 


real 


step length of the electronic iteration 


gamma 


real 


if {i-edyn=2): damping parameter 




real 


secoiiQ step leiigxii oi eiecTjioiiic ixeiaxion cpb-Ciig-aioj 


gam,ma2 


real 


second damping parameter (c.f. gamma and eps-chg-dlt) 


eps-chg-dlt 


real 


if total energy varies less than eps-dlt-chg, delt2 and gamma2 
replace delt and gamma (is reset after moving nuclei) 




i cdi 


liiUlioL. Llldi Li y lidlillL/O . Llllic o LCU iUi tlitJ iiltC^i dUlvJil Ui tliC 

equations of motion (a.u.) 


nOrder 




if (^c/^/n=0,l): order of the scheme: 
predictor corrector 







1 2 




1 


2 3 











3 


4 6 


pfft-store 


real 


percentage of wave functions for which a second transforma- 

UiUli tU ItJdi bJJdCc lo dVUlLlcU. ICi. tt/A.1 I 


iiiooii—Uji^ciii iicy 


1 t?cli 


U-C^icc UU WlilCll Lilt; odlllJJillig, llicUiclil olidil Vjxi odtloilcLl IC.i. 

text) 


idyn 




scheme for solving the equation of motion of the nuclei: 







predictor-corrector 




1 


predictor 
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inp.mod continued 



parameter 


type/range 


explanation 




2 


Verlet 


i pdiiTi 




srheme to itpr?ite the w?ivp fiTnctinns* 







steenest ciescent 




1 


Williams-Soler algorithm 




2 


damped Joaimoupolos algorithm 


% xc 




— Til npfioTiiil * 







LDA (Cerperly/Alder, Perdew/Zunger) 




1 


Becke-Perdew XC fBP) 




2 


Perdew (PW91) 


t-postc 


logical 


.true.: post LDA with functional Lxc 
.false.: start with functional Lxc 


trane 


logical 


.true.: perturbe initial wave functions 


ampre 


real 


if (trane=. true.): amplitude of random perturbation added 
to wave function. 


tranp 


logical 


.true.: perturbe positions of nuclei that are al- 


lowed to move (c.f. tford and amprp) 


avnprp 


real 


if (iranp=. true.): amplitude of random perturbation of 
atomic positions (in bohr) 


tfor 


logical 


.true.: structure optimization (c.f. epsfor, 
force.eps, and, tford), set tdyn^.ialse. 


tdyn 


logical 


.true.: molecular dynamics (c.f. force_eps and 
Iford). sol IJor=. [iihc. 


tsdp 


logical 


structure optimization scheme: 

.true.: modified steepest descent scheme 

.false: damped dynamics scheme 


nstepe 


integer 


if {tfor or tdyn=.tiue.): maximum number of electronic it- 
erations allowed to converge forces, the program terminates 
after nstepe iterations (c.f. force^eps) 


tdipol 


logical 


.true.: employ surface dipole correction. 






the self-consistence cycle is terminated, if for the last three 






iterations 


epsel 


real > 


the variation of the total energy is less than epsel 


epsekinc 


real > 


and, if the average change of wave functions is less than 



epsekinc 
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inp.mod continued 



parameter type/range explanation 

epsfor real > and, if (tfor=.true.), forces on ions with tford^.tiue. are 

smaller than epsfor 



force_eps 2x real > convergence criteria for local and total forces: 

force-eps{l): maximum allowed relative variation in 

local forces before, if {tfor=. true.), exe- 
cuting a geometry optimization step or 
if {tdyn=. true.) calculating total forces. 



force-eps{2): maximum allowed relative variation 
in total forces before moving ions 
(id?/n=.true.) 



max-no-force 


integer 


maximum number of electronic iterations for which no local 
forces shall be calculated per atomic step 


init_basis 




type of basis-set of initialisation: 




1 


plane-wave basis-set (c.f. ecuti) 




2 


LCAO basis-set (c.f. tinit_basis) 




3 


mixed basis-set: LCAO and plane waves 
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Table 3.2. input file start.inp 
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start. inp 



parameter 


type/range 


explanation 


nsp 


integer 


number of atomic species 


neLexc 


real 


number of excess electrons 


nempty 


integer > 


number of empty states 


ibrav 


integer 


cell type (c.f. celldm and latgen.f), e.g.: 




1 


simple cubic lattice, 




2 


fcc-lattice 




3 
8 


bcc-lattice 
orthorhombic 


pgind 


integer 


point group: 







automatic (symmetries and center) 




1 


no symmetries 

subgroup of lattice pomt group 
(c.f. lorav, lib.sym and latgen.f ) 


celldm(1..6) 


6 X real 


lattice parameters of the super cell (depends on ibrav). 
Usually cellam(lj contams the lattice constant m bohr 


nkpt 


integer 


number of k-pomts 


xk(1..3), wkpt 


nkpts lines of 
4 X real 


k-points (c.f. t^relative) and weights of the k-point, the 
sum over which must be 1 


Lfacs(1..3) 


3 X integer >0 


k-point folding factors (c.f. text) - 1 1 1 : no folding 


t-kpointjrel 


logical 


.true.: frame of reference for k-points is 

J "U J-l • 11 J-J-' J- 

spanned by the reciprocal lattice vectors 

.false.: k-points are given in Cartesian coordi- 
nates in units of 2t: jaiat 


ecut 


1 

real 


plane-wave energy cut-on (m Kydj 


ecuti 


real 


plane-wave energy cut-off of the initial wave function 


ekt 


real 


temperature of the artificial Fermi smearing of the elec- 
trons in eV (c.f. tmetal) 


tmetal 


logical 


occuppy electronic states according to: 
.true.: a Fermi distribution (c.f. ek^} 



start. inp continued 



parameter type/range explanation 

.false.: a step-like 

distribution 



tdeqen logical 



.true.: 0(^pation numbers are read from 



start. inp continued 


parameter 


type/range 


explanation 


nfi^rescale 


integer >0 


if {nthm^l): number of time steps before velocities are 
rescaled 


tpsmesh 


logical 


pseudo potentials are: 

.true.: tabulated on logarithmic mesh 

.false.: set up from parametrized form 


coordwave 


logical 


.true.: if {nrho=2), positions of nuclei are read 
from fort. 70 


nsp records describing each atomic species are expected below 


parameter 


type/range 


explanation 


na 


integer 


number of ions of this species 


zv 


integer 


valence charge 


QitOTTl 


cnai dc lei i u 


name oi ine pseucio poieiiLiai 


rgauss 


real 


radius of Gaussian charges 


ionjac 


real>0 


if {tfor=. true.): mass parameter 

if {tdyn=. true.): mass of the nuclei in [amu] 


ion-damp 


l>real>0 


if (i/or^ .true, and tsdp— .false. ): damping parameter 


Lmax 


1,2,3 


highest angular momentum of the pseudopotential (1: s, 
2: p, 3: d) 


IJoc 


integer 
< Ijnax 


angular momentum of the local pseudo potential 


tJnitJ)asis 


3 X logical 


.true.: include s,p and d orbital rsp. in LCAO 
basis-set (c.f. init_basis) 


tauO(1..3), 


3 X real. 


coordinates of the nuclei (units depend on ibrav), flag 


tford 


logical 


whether nuclei may move 


vauO(1..3) 


3x real 


atomic velocities in a.u., expected only if {tdyn=. true. 



and npos—l,A) 
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Table A.l. General symbols and the corresponding variable names. 



symbol 


variable 








number or species 




na(is) 


number or atoms per species 


n 


n,nx 


number ot states 


ns 


nrot 


number ol point symmetry operations 


^loc 


11 1 1 1 

IJOC = /loc + 1 


angular momentum or local pseudopotential 


7 

'max 


7 7 _i_ 1 
l-TflCLX — ^max 


1 • 1 J- 1 J- 

nignest angular momentum 


S 


5(1-3, 1-3, zrorj 


point symmetry operation 


«lat 


alat 


lattice constant 


Ci,k+G 




wave lunction coemcients 




JOCC[l,lK) 


occupation numbers 




wik 


weights of k-points 




Omega 


volume of super cell 




tau0{l-3,ia,is) 


ionic coordinates 
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Table A. 2. Range over which the summation associated with the tabulated index is carried out. 



index 


range 


variable 


summation 


is 


l<is< risp 


is 


species 




-1 ^ • atom 

1 <la< nf^°'^ 


ia 


ions or a species is 


i 


1 < t < n 


i 


electronic states 


I 


^ / ^ ^max) ^ 7^ ^loc 


iJm 


angular momentum 


m 


—l<m<l 




quantum numbers 


G 


1 1 1 1 1 1 1 2 ^ TP 

2 1 + K| 1 S -C/cut 




1 1 J- J-' J- 

reciprocal lattice vectors 


G 


i||G|P<4£;,,t 


ig 




k 


as given in input files 


ik 


k-points 


s 


all point symmetries 


irot 


point symmetry operations 
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Table B.l. parameter file parameter.h 



parameter.h 



parameter 


type/range 


explanation 


nsx 


integer 


maximum number of ionic species 


nax 


integer 


maximum number of ions per species 


nx 


integer 


maximum number of electronic states per k-point 


ngwx 


integer 


maximum number of plane-waves, 


ngw 


'&xngwx 


maximum number of courier components of e.g. the elec- 
tron density 


nxJnit 


integer 


maximum number of basis functions m the initial diag- 
onalization 


ngwix 


integer 


maximum number of plane waves in the initial diagonal- 
ization 


nx-basis 


integer 


maximum number of atomic orbitals in the initial diag- 
onalization 


max-basis-n 


integer 


max{nx, nxJ)asis) 


nlmaxJnit 


integer 


maximum number of atomic orbitals per atom in initial 
diagonalization 


nrlx, nr2x and 
nrSx 


integer 


size of the Fourier mesh (x, y and, z component rsp.) 


nnrx 


(nrix-/--/) xnrSx number of mesh points 




xnrSx 




nkptx 


integer 


highest number of k-points 


ri-fft-store 


1 < integer 


number of wave functions for which a second transfor- 
mation to real space is avoided (c.f. text) 


nlmax 


integer 


nlmax > ElfS"''"^(2 i + 1) - 2lJtoc - 1 


mmax 


integer 


dimension of pseudo potential grid 


nschltz 


integer 









scheme Ledyn^2 disabled 




1 


all schemes enabled 



proper setting reduces storage needs 
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Table B.2. input file inp.ini 
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inp.ini 



parameter 


ty])e/range 


ex])lanation 


ibrav, pgind 


2 X integer 


c.f. start. inp 


nel, tmetal, ekt and, 


real, logical, real 


nel: number of electrons 


tdegen 


and, logical 


c.f. start. inp 


ecut, ecuti 


2 X real 


c.f. start. inp 


tmold, tband, nrho 


2 X logical, integer 


c.f. start. inp 


npos, nthm, nseed 


3 X integer 


c.f. start. inp 


TJon, TJnit, Q, 


3 X real, integer 


c.f. start. inp 


nfi_rescale 






nsp, tpsmesh, 


integer, logical 


c.f. start. inp 


coordwave 






nkpt 


integer 


c.f. start.inp 


xk(l-3), wkpt 


nkpt lines of 






3 X real, integer 


k-pomts m Cartesian co- 
ordinates (units: 27r/aiat) 


al(l-3) 


3 X real 


lattice vectors in a.u. 


a2(l-3) 


3 X real 




aSfl-S) 


3 X real 




bl(l-3) 


3 X real 


reciprocal lattice vectors 


b2(l-3) 


3 X real 


in 27r/ aiat 


1,3(1-3) 


3 X real 




alat, omega 


2 X real 


lattice constant and cell 
volume in a.u. 


nsp records describing each ionic 


species (c.f. start.inp) 


name, na, zv, ionjac 


character* 10, 
2 X integer, real 




ion_damp, rgauss, 


2 X real, 2 x integer 




Lmax, IJoc 






t-initJbasis 






tau0(l-3), 

t_auto_coord(l-3), 

tford 


3 X real, 2 X logical 


c.f. start.inp 
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inp.ini continued 


parameter 


ty])e/range 


ex])lana.tion 






t_ciuto_coord'. c.f. text 


vauO{l-?>) 


3 X real 


c.f. start, inp, expected 
only if (^npos—1,4:^ 


ineq^pos 


3 X integer 




nrot 


integer 


number of point symme- 
tries 


nrot 3x3 matrices representing point group elements 


s(irot) 


3 lines of 
3 X integer 


point group clement pre- 
ceded by a dummy line, 
the frame of reference is 
spanned by the lattice 
vectors al, a2 and a3 
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